Analysis of instabilities and pattern formation in time fractional reaction-diffusion 
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We analyzed conditions for Hopf and Turing instabilities to occur in two-component fractional 
reaction-diffusion systems. We showed that the eigenvalue spectrum and fractional derivative order 
mainly determine the type of instability and the dynamics of the system. The results of the linear 
stability analysis are confirmed by computer simulation of the model with cubic nonlinearity for 
activator variable and linear dependance for the inhibitor one. It is shown that pattern formation 
conditions of instability and transient dynamics are different than for a standard system. As a result, 
more complicated pattern formation dynamics takes place in fractional reaction-diffusion systems. 



I. INTRODUCTION 

In reaction-diffusion systems a stable equilibrium so- 
lution usually changes spontaneously with parameters to 
limit cycle by Hopf bifurcation or stationary dissipative 
structures by Turing bifurcation. As a result, we obtain 
nonlinear dynamics leading to stationary or oscillatory 
structures. When conditions of both instabilities arise, 
we can expect more complex dynamics [lj, [2J, l3( • In case 
of a fractional reaction-diffusion (FRD) system dynamics 
can be much more complex 0, i, i, 0, 0, !, ES M, d • 

Last investigation showed that many complex het- 
erogenous systems are described by differential equations 
with fractional derivatives to represent their anomalous 
behavior [H, El EE El ES El • Therefore, investigation 
of pattern formation in FRD system has both theoretical 
and applied interest. 

In this article the FRD system with cubic nonlinearity 
is studied for the case when both instabilities take place. 
We have focused on the dynamics of FRD model under 
conditions when sufficiently complex patterns arise in the 
system dynamics. The results of analytical treatment of 
the linearized model arc validated by computer simula- 
tions of nonlinear dynamics. 



II. MATHEMATICAL MODEL 



Let us consider the FRD system 

T c U t = lu xx + W{u, A), 



(1) 



with two variables u = (ui.u 2 ) T on the x £ (0,£) sub- 
ject to Neumann: u x \ x =o,Cx ~ boundary conditions 
and with certain initial conditions, W — (Wi,W2) T , 
W\,W% - smooth reaction kinetics functions, A - real 
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parameter, r and I are positive diagonal matrices t = 
diag[ri], I — diag[lf\ > . 

Fractional derivatives c u t on the left hand side of the 
equations ([1]), instead of the standard time derivatives, 
are the Caputo fractional derivatives in time [HI, Ell of 
the order < a < 2 and are represented as 



d>(t) 



r(m - 



Am) 



(r) 



< D J (t- 
o 



dr, (2) 



where 



1 < a < m, m € 1, 2. 



III. LINEAR STABILITY ANALYSIS 

Due to the property of Caputo derivative the sta- 
bility of the steady-state solutions of the system 
(HJ) corresponding to homogeneous equilibrium state 
W(uo,Ao) = can be analyzed by linearization 
of the system nearby this constant solution uq — 
(ui,U2) T ■ The linearization of FRD system ([1]) 
leads to fractional ODEs with right hand side matrix 
(an - k 2 h\)/Ti a 12 /n 

a2i/T 2 (a 22 - k 2 lj)/T 2 
form of which is given by eigenvalues Xi 2 = ^(trF ± 

\Jtr 2 F — 4det F) (coefficients atj represent Jacoby ma- 
trix) . 

For a : < a < 2 for every point inside the parabola 
det F = tr 2 F/4, we can introduce a marginal value a : 
ao = —\Arg(Xi)\ given by the formula [1,0] 



F(k) 



diagonal 



a 



| arctan ^4det F/tr 2 F - 1, trF > 0, 
2 - | arctan y/ 4 det F / tr 2 F - 1, trF < 0. 



(3) 



The value of a is a certain additional bifurcation param- 
eter which switches the stable and unstable states of the 
system. At lower a : a < cvo = § the system 
has oscillatory modes, but they are stable. Increasing 
the value of a > ao = ^|^4rf;(Ai)| leads to oscillatory 
instability. 
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It is widely known for integer time derivatives [J 0, 
3] that system (JTJ) becomes unstable according to either 
Hopf (fc = 0) 

trF>0, detF(0)>0. (4) 
or Turing (fco 7^ 0) bifurcations 

trF<0, detF(0)>0, det F(k ) < 0. (5) 

and these both types of instabilities are realized for pos- 
itive feedback (an > 0) [EH 9. 

In the case of fractional derivative index, Hopf bifurca- 
tion is not connected with the condition an > and can 
hold at a certain value of a when the fractional derivative 
index is sufficiently large [l(|. Moreover, in fractional RD 
systems at a > 1 when it is easier to satisfy conditions of 
Hopf bifurcation, we meet a new type of instability [H, Q 

trF<0, 4detF(0) <tr 2 F(0), 4 det F(k ) > tr 2 F(k ) 

(6) 

It is worth to analyze inequalities ([6]) in detail. Taking 
into account explicit form of F(k) the last two conditions 
can be rewritten as: 



(aim - a 2 2T 2 ) 2 > ™4ai2a 2 iTiT2, 



(7) 



-4ai 2 a 2 iTiT2 > [(an - fc 2 Z 2 )r 2 - (a 22 - fc 2 Z 2 )Ti] 2 . (8) 

The simplest way to satisfy the last condition is to esti- 
mate the optimal value of k = fco 



fc = 2 



'012021 



11/T2 - q/n 



1/2 



(9) 



Having obtained ([5]), we can estimate the marginal value 
of «o 



ao = 2 arctan T, 



where the expression T has the following view 



T 



(-4ai 2 a2iTiT2 



,1/2 



(a n T 2 - a 22 Ti 



(10) 



(11) 



anT 2 — a 2 2Ti 



The analysis of expressions (JSJ) shows that at k = we 
have two real eigenvalues that are less than zero, and the 
system is certainly stable for the Hopf bifurcation. If the 
last inequality takes place for a certain value of fco 7^ 0, 
we can get two complex eigenvalues, and a new type of 
instability, connected with the interplay between the de- 
terminant and the trace of F(k) of the linearized system, 
emerges. With such type of eigenvalues, it is possible to 
determine the value of fractional derivative index when 
the system becomes unstable for Hopf bifurcation with 
this wave number [8J]. 



IV. FRACTIONAL REACTION DIFFUSION 
SYSTEM WITH CUBIC NONLINEARITY 



To demonstrate the properties of FRD system, let us 
consider the model with cubic dependance for activator 
variable W% = ui—u\ — u 2 and the linear for the inhibitor 
variable W 2 = —u 2 + f3ui + A. This model was proposed 
firstly by R. FitzHugh |20| for description of the prop- 
agation of voltage impulse through a nerve axon and is 
known as Bonhocffer-van der Pol model. In RD systems 
this model was considered in many books and articles 
(see for example [3, H, Hj] ) . The homogeneous solution of 
variables u± and u 2 can be obtained from the system of 
equations Wi(ui,u 2 ) = 0, W 2 (Tii,u 2 ) = 0, which in their 
turn determine two nullclines u 2 (u\ ) . The intersection of 
these nullclines in the point P = (u\,u 2 ) is determined 



by equation: u 



fiu\ — A =0. In this case the val- 



ues of external parameters A,/3 determine the value of ui 
and this makes it possible to investigate the conditions 
of different types of instability explicitly considering pa- 
rameter u\ as the main parameter for system analysis. 

For investigation of the Hopf bifurcation let us consider 
homogeneous perturbation with k = 0. The linear anal- 
ysis of the system with a = 1 shows that, if ti/t 2 > 1, 
the solution corresponding to any intersections of two 
isoclines is stable. The smaller is the ratio of ti/t 2 , the 
wider is the instability region (solid lines with points on 
Fig. 1). Formally, at t\/t 2 — > 0, the instability region for 
Iti coincides with the interval (—1, 1) where the null cline 
W(ui,u 2 ) — has its increasing part (Fig. 2b). These 
results are very widely known in the theory of nonlinear 
dynamical systems [if H, H[ . 

In the fractional differential equations the conditions 
of the instability depend on the value of a and we have to 
analyze the real and the imaginary parts of the existing 
complex eigenvalues, especially the equation: 

***F-*>F = ^-D+^-f^ _ IV > 0. 

T\T 2 \ Tl T 2 J 

In fact, with complex eigenvalues, it is possible to find 
the corresponding value of a where the condition a > aa 
is true. Omitting simple calculation, we can write an 
equation for marginal values of u\ 



n?-2(l + -)u 2 + ^ -2^(2/3- 1) + 1 = 0, (13) 
t 2 "nf t 2 



and solution of this biquadratic equation gives us the 
domain where the oscillatory instability can arise 



1 



n 

t 2 



± 2, 



(14) 



This expression estimates the maximum and minimum 
values of Mi where the system can be unstable at marginal 
value of a = «o = 2 as a function of t\/t 2 and f3. 

The typical stability domains for considered FRD sys- 
tem in the coordinates (ui,ri/r 2 ) for different values of 
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fractional derivative index a are presented in Fig.la-d, 
where curves corresponding to a = 1 are denoted by 
more solid lines. This makes it possible to see how other 
curves a ^ 1 are located with respect to standard system 
a = 1. The solid points on both sides denote the inter- 
val of maximum instability. For each particular value a 
in the region between the corresponding curve and hor- 
izontal axis, the system is unstable with wave numbers 
k — 0, and outside it is stable. Figures (a) and (b) present 
the plots for the value of (3 = 1.01. The left-hand side 
plot corresponds to a < 1 and the right-hand side plot 
corresponds to a > 1. It is easy to see from the plot 
(a) that if the value of T1/T2 increases (we move along 
vertical axis) the instability domain decreases (a < 1). 
At T1/V2 = 1 and a = 1 it vanishes completely. The 
solid point in the middle of each plots corresponds to the 
case when the system becomes stable for all values of 
u\ : point (0,1) in coordinates (u\,ti/t2). The situation 
changes for a > 1. The system is unstable not only for 
T1/T2 < 1 but also for t±/t2 > 1. An increase in a makes 
the instability domain much wider with respect to two 
coordinates (mi,t 1 /t 2 ) and we obtain butterfly like do- 
mains for a > 1. This means that with increasing t\/t% 
(moving along vertical axis) the system becomes stable in 
the center and unstable at greater values u\ In such case 
the instability domain becomes symmetric along vertical 
axis with the minimum point at the u\ = 0. We observe 
similar behavior for large values of (3. Instability domains 
in Fig. lc,d are presented for j3 = 10 and show the same 
trend with respect to a, but the region for a > 1 is much 
greater than for (3 = 1.01. At the same time, for a < 1 
and (3 = 10 instability domain shrinks very sharply in 
comparison to the same domain plot when (3 = 1.01. 

It is possible to obtain solid understanding of the mech- 
anism of the instability from the plot of eigenvalues. Typ- 
ical instability domains for the same parameters as on 
Fig. la,b for k = are presented on the Fig. 2a. Hori- 
zontal lines (i, ii, iii) on the instability domain plot cor- 
respond to eigenvalues plots below. Let us analyze each 
of the possible situations in more detail. 

For the case (i) we have sub-domains with real posi- 
tive, real negative and complex eigenvalues. The easiest 
way of obtaining instability is realized at |wi | < uf when 
all the roots are real and positive (Fig. 2a(i)). This 
region is presented by dark grey color and positive eigen- 
values mean that the system is unstable practically for 
any value of a > 0. Inside the domain \uf\ < \ui\ < \u^\ 
there is a certain domain of a : (0 < a < 2) where the 
Hopf bifurcation takes place. Point D divides the region 
into two domains where Re A < and ReA > 0. In the 
domain ReA < the system could be unstable according 
to greater values of a > 1. In turn, for ReA > 0, the 
system could be stable at a < 1. In other words, between 
points C and E we have eigenvalues with imaginary part, 
and the value of a can change the stability of the FRD 
system. In the domain \u\ | > |uj | we have two real and 
negative roots and, as a result, the system is stable. 

For system parameters corresponding to the case (ii) 




-1.0 -0.5 0.0 0.5 1.0 -2.0-1.0 0.0 1.0 2.0 




-1.0 -0.5 0.0 0.5 1.0 -8.0 -4.0 0.0 4.0 8.0 



(c) (d) 



FIG. 1: Instability domains in coordinates (ui,ti/t2) for a 
fractional order reaction-diffusion system with sources Wi = 
u\ — uf — U2, W2 = — U2 + f3ui + A for different values of 
ao = 0.25,0.5,0.75,1.0,1.25,1.5,1.75. The results of com- 
puter simulation obtained at h = I2 — for: /3 — 1.01 - (a), 
(b) and j3 = 10.0 - (c), (d) . 



the real part of eigenvalues becomes less than zero for 
all Ui . At the same time, for |ui | < |uf \ the roots are 
complex and and according to condition ([3]) instability 
takes place for a > ao > 1. For |7ii| > \u^\ the roots 
become real and negative, and the system is stable. 

In the case (iii) at the center of u\ ( \u\\ < \uf \ ) we 
have two real negative eigenvalues, and the system is 
stable. For |7i^| < |7ii| < \u^\ we have complex roots 
and certainly according to condition instability takes 
place for a > olq > 1. In this case, the instability domain 
consists of two symmetrical regions separated by a stable 
region at the center where the system is stable for any a. 
For \u\ I > Hi the system is stable again. 

Let us analyze the Turing Bifurcation {k ^ 0). Eigen- 
values for different values of k are presented in Fig. 2b. 
The top plot corresponds to nullclincs of the systems just 
to show that nullcline intersection determines eigenval- 
ues in Fig. 2b(iv). We present eigenvalues for k = 1 
and k = 2 and for comparison k — 0. It can be seen 
from the picture (iv) that at intersection of nullclines in 
the vicinity of zero value of Hi nonhomogeneous modes 
have much greater values and we can expect a forma- 
tion of stationary dissipative structures. If the ratio I1/I2 
is sufficiently small, Turing bifurcation is dominant for 
all region Hi < 1. Analyzing (O we can conclude that 
these conditions are practically the same for fractional 
and standard RD systems. However, what is very im- 



FIG. 2: Instability domains (ReX - black lines, ImX - grey 
lines) for k = 0, j3 = 1.05 and different proportions of T1/T2 = 
0.5- (z),1.0- (u),2.0- (Hi). 

portant is that the transient processes and the dynamics 
of these systems are different, and for this reason final 
attractors can often be different even though the linear 
conditions of instability look the same. 

Now let us consider that the system parameters are 
close to the ones represented by point P in Fig. 2b. As 
a result, for certain ratio of ImX and ReX we expect the 
formation of oscillatory inhomogeneous structures. How- 
ever, if the solution ui is close to zero, the decrease of a 
leads to steady state dissipative structures. Such trend 
is quite general and if in standard system we have steady 
state solutions, the increase of a in FRD system leads 
to non stationary structures. In this case, by changing 
intersection point of nullclines or value of a we can stim- 
ulate stationary or temporary pattern formation. If the 
absolute value of eigenvalues for k — and k ^ are 
comparable we can expect more complex spatio-temporal 
dynamics. 

Above we have considered that the linearized system 
is unstable for either Hopf or Turing bifurcation. Below, 
we consider the case when we don't have Turing or Hopf 
bifurcation. For realization of instability conditions ([5]) 



FIG. 3: The null-clines for /3 — 2.1, A — 0.5 and eigenvalues 
for different values of k (k = - hair-lines, k = 1 - dash lines, 
k — 2 - thick lines) -(b). The eigenvalues are presented for 
the following parameters: l\ — 0.025, /3 = 2.1,ti/t2 = 0.21 
- (iv), h = 0.1, P = 1.01, n/ra = 0.6 - (v), h = 2.1, (5 = 
1.01, n/r 2 = 3.5 - (vi). 



for k 7^ 0, ImX ^ the fractional derivative index must 
be greater than some critical value ao- Eigenvalues for 
such instability are presented in Fig. 2b(v,vi). We can 
see that outside a small domain in the center the system 
is stable for k — 0. At the same time on this interval 
we have complex eigenvalues for k ^ (Fig. 2b (v)). In 
the plot (v) for l\/l2 < 1 we have a separate domain for 
k = 2 where we can expect inhomogeneous oscillations 
with this wave number. 

In Figure (vi) we present the situation where all roots 
have ReX < and in standard system we do not have 
instability at any values of tlx . In FRD system the roots 
are complex for select values of k (k ^ 0). In other words, 
for such system at a > 1 we can obtain conditions of Hopf 
bifurcation (O which lead to inhomogeneous oscillatory 
structures [8|, [9( even for I1/I2 > 2. This situation can be 
predicted from symmetrical view of expression (jTTJ) for 
the system under consideration 
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(a) 



(b) 





(c) 



(d) 



FIG. 4: The view of the surface T in coordinates 
(h/h, T2/T1) for {3 — 2 and different values of u\ (iti =0.1 
- (a), iti = 5.0 - (b), ui = 1.25 - (c), ui = 1.5 - (d)) 



(1 



«fr 2 + Z|n 



1(t2 - l%Ti 



(1 



r 2 /ri + 1 



(15) 

The plot of these surfaces, as a function of h/h and 75/71, 
are presented in plots 3 (a-d) for different values of u\ . 
We can see that in Fig. (a) and (b) the maximum value 
of T is reached at the boundary and h/h « 1 (Fig. 3a) 
and t 2 /ti « 1 (Fig. 3b). For figure. 3 c,d the optimal 
instability conditions are reached at certain combination 
of the parameters h/h and T1/T2 and we can expect inho- 
mogeneous oscillations at different relationships of I1/I2 
and ri/r 2 . 

We expect the oscillatory structures to emerge if at 
given values of Hi, h/h and ti/tz, the fractional deriva- 
tive index a is greater than the one represented on the 
surface and less than the one needed for Hopf bifurca- 
tions for k = 0. This means that only the perturbations 
with these wave numbers are unstable, and they are un- 
stable for oscillatory fluctuations. This situation is qual- 
itatively different from the integer RD system whether 
either Turing (k 7^ 0) or Hopf bifurcation (k = 0) takes 
place, and this depends on which condition is easier to 
realize. Thus, in the system under consideration, we can 
choose the parameter when we don't have Turing and 
Hopf bifurcations (for k = 0) at all. Nevertheless, we ob- 
tain that conditions for Hopf bifurcation can be realized 
for a nonhomogeneous wave number. As it is seen from 
the figure, there are conditions where only instability ac- 
cording to non-homogeneous wave numbers holds. As a 
result, perturbations with k = relax to the homogenous 
state, and only the perturbations with a certain value of k 



become unstable and the system exhibits inhomogeneous 
oscillations. 



V. PATTERN FORMATION 

The results of the numerical simulation of the frac- 
tional RDS (P} are presented on Fig. 4,5. From the 
pictures we can see that in such system we obtain a rich 
scenario of pattern formation: standard homogeneous os- 
cillations, Turing stable structures, interacting inhomo- 
geneous structures and inhomogeneous oscillatory struc- 
tures. We have obtained that the ratio of characteristic 
times and the order of fractional derivative qualitative 
transform pattern formation dynamics: homogeneous os- 
cillations in the first limiting case and stationary dissipa- 
tive structures in the second one [l(| ■ Here we show that 
the change in any parameter which qualitatively changes 
the eigenvalues of the linearized system can change the 
system dynamics. Spatiotemporal dynamics of the FRD 
system can mainly be determined by the maximum eigen- 
values for the corresponding modes. In Fig. 4a, b we can 
see stationary dissipative structures as a result of for- 
mation of the unstable mode presented in Fig. 2b (iv) 
for a < 1. External parameter A determines the inter- 
section point and the slope of the isoclines in this point 
and the power for each particular mode at this parame- 
ter. In particular, for A — 0.25 nullclines intersect at the 
point where maximum value has eigenvalue with k = 2, 
which is responsible for Turing bifurcation with station- 
ary structures. For A=0.5 (Fig. 2b(iv)) nullclines inter- 
sect at the point where Hopf bifurcation takes place. The 
characteristic feature for these two limit cases is the in- 
stantaneous formation of either dissipative structures or 
homogeneous oscillations. Increasing influence of Hopf 
bifurcation when the Turing one is dominant, or increas- 
ing Turing bifurcation when Hopf is dominant, leads to 
more complicated transient dynamics (Fig. 4c, d). When 
conditions of these two instabilities practically coincide, 
we can obtain either oscillatory inhomogeneous struc- 
tures or modulated homogeneous oscillations (Fig. 4e,f). 
Moreover, at parameters, when the real part of eigen- 
values is close to zero, small variation of a changes the 
type of bifurcation. This trend is typical for any a < 1. 
For a > 1 the structure formation can be much more 
complicated. 

Let us consider the bifurcation diagram presented in 
Fig. lb,d. It was already noted that the region inside 
the curve is unstable for wave numbers k = and out- 
side - it is stable. From the viewpoint of homogeneous 
oscillations, the system is stable near u\ — 0. However, if 
we have l\ << I2, the system becomes unstable according 
to Turing instability. As a result, we expect the forma- 
tion of stationary inhomogeneous structures. In fact, at 
the beginning, only inhomogeneous fluctuations grow in 
amplitude and lead to inhomogeneous pattern formation. 
At the same time, at the dynamics of structure forma- 
tion, the amplitude of the structures increases, and at 
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FIG. 5: Dynamics of pattern formation for u\ variable. The 
results of computer simulations of the system at parameters: 
a = 0.8, A = -0.25, P = 2.1, h = 0.025, h = 1, ti/t 2 = 0.1 
- (a); a = 0.8, A = -0.55, P = 2.1, h = 0.025, h = 1, 
ri/r 2 = 0.1 - (b); a = 0.8, A = -0.4, p = 2.1, h = 0.025, 
h = 1, n/r 2 = 0.1 - (c); a = 0.8, A = -0.45, /3 = 2.1, 
/i = 0.025, h = 1, ri/r 2 = 0.1 - (d); a = 1.6, A = -0.01, 
P = 1.05, «i = 0.05, h = 1, Ti/r 2 = 1.45 - (e); a = 0.7, 
A = -0.3, P = 2.1, /i = 0.05, h = 1, n/r 2 = 0.2 - (f); 



maximum and minimum amplitude, the structures fall 
into the domain where the homogenous structures are 
unstable. As result, we obtain complex interaction of 
Turing and Hopf bifurcations (Fig. 5a, b). 

Inhomogeneous oscillatory structures are presented on 
pictures (Fig. 5 c,d)). Due to different evolution pattern 
for two variables, the activator variable u\ is presented 
in the left column and the inhibitor one U2 is presented 
on the right column. Such structures are obtained at 
eigenvalues with negative real part when the intersection 
of the null-cline for activator variable is located on the 
decreasing part of null cline U2{u\). The corresponding 
eigenvalues for this situation are presented in Fig. 2b (v). 
From the eigenvalue plots (Fig. 2) we can see correspond- 
ing separated domains where inhomogeneous oscillatory 
modes with k = 2 are unstable. Successive increase of 
the fractional derivative index will increase the amplitude 
of the presented in Fig. 5a, b inhomogeneous oscillation. 
As a result, ingomogeneuos oscillatory structures of large 



FIG. 6: Dynamics of pattern formation for ui (left column) 
and M 2 (right column) variables. The results of computer 
simulations of the systems at parameters: A = —0.01, a = 
1.8, P = 1.01, h = 0.02, h = 1, ri/r 2 = 3.5 - (a-b); A = 1.95, 
a = 1.82, = 1.01, h = 0.1, h = 1, ti/t 2 = 0.6 - (c-d); 
A = -0.01, a = 1.75, P = 10, h = 0.05, h = 1, n/r 2 = 0.05 

- (e-f); 

amplitude are realized in the system. 

Another type of oscillatory structures with a little bit 
more complicated dynamics are presented in Fig. 5e,f for 
P = 10, T\ « t 2 and a ;$ 2. As in the case considered 
above we have inhomogeneous oscillatory structures the 
surface of which in the region of slow motion oscillates 
with fast frequency 1 jr\ . Such behaviors are due to oscil- 
latory property of the activator system at a approaching 
the value of 2. 



VI. CONCLUSION 

We have shown a complex spatio-temporal pattern for- 
mation in simple FRD system and compared these results 
with standard one. The fractional derivative index plays 
a crucial role in this pattern formation because conditions 
of time bifurcations depend substantially on its value. By 
eigenvalue analysis we have studied instability conditions 
for k = 0, and k ^ for different values of external pa- 
rameter A (the same as u\). Nonlinear solutions show 
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that dynamics of the system is determined by most un- 
stable modes. When linear increments are comparable, 
we have an interplay between Hopf and Turing modes 
leading to more complex dynamics. 

When a fractional derivative index is changed from to 
1, the large-amplitude structures are stationary if a limit 
cycle is damped at n ~ r 2 . Oscillatory structures at 
these values of fractional derivative index can be realized 
only at ri << t 2 . 



When a fractional derivative index is changed from 
1 to 2, the large- amplitude structures have more com- 
plex dynamics. Moreover, the spatiotemporal structures 
are observed even at n >> r 2 . Complex structures are 
observed in the region, when the bifurcation parameter 
leads to Turing and Hopf instabilities, as well as in the 
regions where these instabilities are damped. The system 
moves to large amplitude limit cycle with further change 
in the fractional derivative index. 
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